function [t,Tg,Tgb,ATg,N,Dcum,Nd,Ns,Ndead,r,t_b_trat]=kinetic_model(param,dose,halflife,number_of_injections)        

%This code implements a thyroid response model with different types of cells: undamaged, sublethally damaged, lethally damaged and dead, 
%as well as two tumor markers, namely the serum concentration of Tg and TgAb. A pretreatment phase simulates the growth of tissue. This allows to
%obtain a thyroid tissue with a mass equal to some predetermined value (for example the mass of a experimentally measured thyroid/ lesion)
%and the corresponding concentration of serum Tg and TgAb. The treatment phase simulates a I-131 irradiation producing sublethal and lethal 
%damage. The treatment can involve several I-131 administrations, delivered at periodic time intervals. 

%Input: 
%1-param: vector with the model parameter values (see lines 28-48 for a description).
%2-dose: vector with the dose achieved in each treatment administration. 
%3-halflife: vector with the effective half-life of the radionuclide in each treatment adminstration.
%4-number_of_injections: number of treatment administrations.

%Output:  
%1-t:        time vector
%2-Tg:       vector with the temporal evolution of free Tg serum concentration in ng/ml
%3-Tgb:      vector with the temporal evolution of Tg serum  bound to TgAb concentration in ng/ml
%4-ATg:      vector with the temporal evolution of the TgAb serum concentration in ng/ml
%5-N:        vector with the temporal evolution of the number of undamaged cells
%6-Dcum:     vector with the temporal evolution of the absorbed dose 
%7-Nd:       vector with the temporal evolution of the number of doomed cells
%8-Ns:       vector with the temporal evolution of the number of sublethally damaged cells
%9-Ndead:    vector with the temporal evolution of the number of dead cells
%10-r:       vector with the temporal evolution of the dose rate
%11-t_b_trat: time before treatemtn

%Description of input parameters vector
%Proliferation parameters:
lambda=log(2)/param(1);     %Normal cells prolif rate (day^-1)
lambda_d=log(2)/param(2);   %Doomed cells prolif rate (day^-1)
lambda_s=log(2)/param(3);   %Sublethally damaged prolif rate (day^-1)
lambda_t=param(4);          %Thyroglobulin production rate in ng/(ml·day)
%Radiosensitivity, repair & resorption rates
a=param(5);                   %Letal damage (damaged cells/decay)
b=param(6);                   %Subletal damage (damaged cells/decay)
mu=log(2)/(param(7)/24);      %Subletal damage repair rate (day^-1),(param(7) is in hours, 2h halflife)
gamma=log(2)/param(8);        %Resorption rate for doomed cells (Ln(2)/T_1/2, day^-1)
ke=log(2)/param(9);           %Rate of thyroglobulin elimination (day^-1)
N0=param(10);                 %Number of undamaged cells just before treatment
cte_Tg=param(11);             %Thyroglobuline production modulation constant (ng/ml)
Nmax=param(12);               %Threshold for saturation of cell proliferation
gamma2=log(2)/param(13);      %Rate of death for doomed cells (Ln(2)/T_1/2, day^-1)
k=log(2)/param(14);           %TgAb natural elimination rate (day^-1)
ka=log(2)/param(15);          %Free Tg binding rate to Tg antibodies in ml/(IU·day), or ml/(ng day) for Tg antibodies
z=log(2)/param(16);           %TgAb production rate IU/(ng·day)
keb=log(2)/param(17);         %Bound Tg enhanced metabolic clearance rate (day^-1)
zb=param(18);                 %TgAb production rate (lymphocyte related) in IU/(ml·day)
zt=param(19);                 %TgAb production rate in ng/(ml·day)

%Initial Conditions: the simulation begins with a small number of undamaged proliferating cells and all other variables =0

N(1)=100000;     %Initial undamaged cells
Nd(1)=0;         %Initial doomed cells
Ndead(1)=0;      %Initial dead cells
Ns(1)=0;         %Initial subletally damaged cells
Tg(1)=0;         %Initial Tg (free compartment)
Tgb(1)=0;        %Initial Tg (bound compartment)
ATg(1)=0;        %Initial TgAb
r(1)=0;          %Initial dose rate
Dcum(1)=0;       %Initial absorbed dose
t(1)=0;          %Initial time


dt=20/(60*24);    %Time step, 20 minutes in days
resting_time=180; %Time between treatment adminstrations injections (in days)
i=1;
decay_rate=0;    %This is a dummy value for the growing phase. The parameter value will be reasigned in the treatment phase (line 144).


%Simulation of tissue evolution until the thyroid mass equal to the mass achieved when the treatment begins.
while N(end)<N0
      %Normal (undamaged) cells compartment
      dN  = lambda*N(i)*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax)...
           +mu*Ns(i)...
           -a*r(i)*N(i)...
           -b*r(i)*N(i);
       
     %Doomed cells compartment 
     dNd = lambda_d*Nd(i)*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax)...
           +a*r(i)*N(i)...
           +a*r(i)*Ns(i)...
           +b*r(i)*Ns(i)...
           -gamma*Nd(i);
       
     %Dead cells compartment
     dNdead = gamma*Nd(i)...
           -gamma2*Ndead(i);
       
     %Subletally damaged cells compartment 
     dNs = lambda_s*Ns(i)*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax)...
           +b*r(i)*N(i)...
           -mu*Ns(i)...
           -a*r(i)*Ns(i)...
           -b*r(i)*Ns(i);
       
     %Thyroglobulin (free) compartment
     dTg =lambda_t*(N(i)+Ns(i)+Nd(i))/(cte_Tg+(Tg(i)+Tgb(i)))...
          +zt*Nd(i)...
          -ke*Tg(i)...
          -ka*ATg(i)*Tg(i);
      
     %Thyroglobulin (bound to TgAb) compartment
     dTgb=ka*ATg(i)*Tg(i)...
         -(ke+keb)*Tgb(i);
     
     %AntiThyroglobulin autoantibodies compartment
     dATg=z*(Tg(i)+Tgb(i))...
         +zb*Nd(i)...
         -ka*ATg(i)*Tg(i)...
         -k*ATg(i);
     
    %dose rate (Equal to zero. It is defined in this growing phase to have all variable vectors with same size)
     dr  = -decay_rate*r(i);
     
     %Variables update
     N(i+1,1)=N(i)+dN*dt;
     Nd(i+1,1)=Nd(i)+dNd*dt;
     Ndead(i+1,1)=Ndead(i)+dNdead*dt;
     Ns(i+1,1)=Ns(i)+dNs*dt;
     Tg(i+1,1)=Tg(i)+dTg*dt;
     Tgb(i+1,1)=Tgb(i)+dTgb*dt;
     ATg(i+1,1)=ATg(i)+dATg*dt;
     r(i+1,1)=r(i)+dr*dt;
     Dcum(i+1)=Dcum(i)+r(i)*dt;
     t(i+1)=t(i)+dt;
     i=i+1;
end
%Construction of the a vector with the times of treatment administration.
t_b_trat=t(end); %time before treatment
t_end=t_b_trat+(number_of_injections-1)*resting_time+4*30; %after de last treatment adminsitration, the evolution the tissue is simulated for 4 months 
treatment_time=[0 [t_b_trat:180:t_b_trat+180*(number_of_injections-1)] t_end]; 
%Variable used to control the number of treatment administrations already delivered
injection=1;

%Simulation of the treatment
for j=2:(length(treatment_time)-1)
    t_span=(treatment_time(j):dt:treatment_time(j+1)-dt);
    t=[t t_span];
    steps=length(t_span);

    if injection<=number_of_injections
        %Initialization of the dose rate value at the beginning of the adminstration
        decay_rate=log(2)/halflife(injection);
        r0=dose(injection)*decay_rate;
        r(end,1)=r(end,1)+r0;
    end
%Simulation of the temporal evolution of the thyroid and dose variables in the time between administrations
    for i=length(N):length(N)+(steps)-1 
      
        %Normal (undamaged) cells compartment
        dN  = lambda*N(i)*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax)...
           +mu*Ns(i)...
           -a*r(i)*N(i)...
           -b*r(i)*N(i);
 
        %Doomed cells compartment
        dNd = lambda_d*Nd(i)*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax)...
           +a*r(i)*N(i)...
           +a*r(i)*Ns(i)...
           +b*r(i)*Ns(i)...
           -gamma*Nd(i);
       
        %Dead cells compartment 
        dNdead = gamma*Nd(i)...
           -gamma2*Ndead(i);

        %Sublethally damaged cells compartment
        dNs = lambda_s*Ns(i)*(1-(N(i)+Nd(i)+Ns(i)+Ndead(i))/Nmax)...
           +b*r(i)*N(i)...
           -mu*Ns(i)...
           -a*r(i)*Ns(i)...
           -b*r(i)*Ns(i);
       
        %Thyroglobulin (free) compartment     
        dTg = lambda_t*(N(i)+Ns(i)+Nd(i))/(cte_Tg+(Tg(i)+Tgb(i)))...
          +zt*Nd(i)...
          -ke*Tg(i)...
          -ka*ATg(i)*Tg(i);
        %Thyroglobulin (bound to TgAb) compartment
        dTgb=ka*ATg(i)*Tg(i)...
         -(ke+keb)*Tgb(i);

        %AntiThyroglobulin autoantibodies compartment
        dATg=z*(Tg(i)+Tgb(i))...
         +zb*Nd(i)...
         -ka*ATg(i)*Tg(i)...
         -k*ATg(i);%-k2*ATg(i)*Tg(i)
        %dose rate
        dr  = -decay_rate*r(i);%*volume_correction;
     
        %Variables update    
        N(i+1,1)=N(i)+dN*dt;     
        Nd(i+1,1)=Nd(i)+dNd*dt;     
        Ndead(i+1,1)=Ndead(i)+dNdead*dt;     
        Ns(i+1,1)=Ns(i)+dNs*dt;     
        Tg(i+1,1)=Tg(i)+dTg*dt;     
        Tgb(i+1,1)=Tgb(i)+dTgb*dt;     
        ATg(i+1,1)=ATg(i)+dATg*dt;     
        r(i+1,1)=r(i)+dr*dt;     
        Dcum(i+1)=Dcum(i)+r(i)*dt;
 
    end
    %Update the administration control variable
    injection=injection+1;
end
end